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LIMITS  OF  APPLICABILITY 
OF  THE  BARBER  AND  HILL  T-Matrix  CODE 


1 .  INTRODUCTION 

A  recent  book  by  Barber  and  Hill*  provides  a  thorough 
and  rigorous  mathematical  treatment  of  light  scattering  from  four 
different  types  of  systems:  plane  slabs,  infinite  circular 
cylinders  illuminated  at  normal  incidence;  axisymmetric 
particles;  and  spheres.  Computer  diskettes,  formatted  for  DOS 
machines  and  packaged  with  the  book,  contain  extensively 
documented  FORTRAN  source  codes  for  the  numerous  programs 
expounded  within  the  text.  In  the  field  of  light  scattering,  it 
is  rare  to  find  computer  programs  in  such  a  polished  and  final 
state;  that  combined  with  the  author's  reputations  for  meticulous 
work  should  insure  that  these  programs  remain  standard  tools  for 
some  time  to  come. 

This  report  is  concerned  only  with  the  chapter  on 
axisymmetric  particles,  which  deals  with  light  scattering  from 
spheroidal  bodies,  both  prolate  and  oblate.  The  Barber  and  Hill 
programs  will  not  work  for  spheroids  that  are  too  large  or  too 
eccentric;  but,  the  text  does  not  adequately  indicate  where  the 
limits  lie.  The  object  of  this  report  is  to  show  the  range  of 
particles  that  can  be  handled  by  the  Barber  and  Hill  routines, 
and  to  suggest  suitable  values  for  the  convergence  parameters 
discussed  below,  which  the  user  must  supply  to  the  computer 
program  for  each  particle. 


2.  OPERATION  OF  T1.F0R 

Nine  programs,  named  Tl.FOR  through  T9. FOR,  are 
supplied  for  calculating  scattering  properties  of  spheroids.  Of 
those,  Tl.FOR  is  the  primary  program  in  that  it  calculates  a 
particle's  T-matrix,  the  mathematical  structure  that  contains 
information  about  the  particles 's  scattering  behavior.  The  T- 
matrix  is  a  function  of  a  particle's  size,  aspect  ratio,  and 
refractive  index,  but  does  not  depend  on  the  particle's 
orientation.  It  is  used  by  the  other  programs  to  calculate 
angular  scattering  intensities  and  absorption,  extinction,  and 
scattering  cross  sections  under  various  circumstances.  (Program 
Tl.FOR  can  also  be  used  to  calculate  a  related  A-matrix,  required 
for  internal  field  calculations  performed  by  program  T8.F0R.) 


♦Barber,  P.W.,  and  Hill,  S.C.,  Light  Scattering  bv  Particles: 
Computational.. Methods ,  World  Scientific  Publishing  Co.  PTE.  LTD, 
Singapore ,  1990 . 
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The  operator  must  run  T1  at  least  three  separate  times 
to  determine  the  convergence  parameters  nrank,  ntheta,  and  nm, 
and  to  produce  the  T-matrix.  Nrank,  which  may  not  exceed  25  in 
this  implementation,  is  the  order  of  the  highest  term  to  be  kept 
in  the  vector  spherical  harmonic  expansion  of  the  electric 
fields.  Nm  is  the  number  of  azimuthal  modes  required  in  the 
expansion,  and  ntheta  is  the  number  of  sample  points  required  for 
accurate  numerical  integration.  As  T1  is  presently  written, 
ntheta  may  not  exceed  100.  Each  time  T1  is  run,  the  user  is 
asked  to  enter  particle  properties  (size  parameter,  axial  ratio, 
and  complex  refractive  index),  values  for  nrank  and  ntheta,  and  a 
case  indicator  (ic)  that  "tells"  the  program  what  test  to  carry 
out. 


In  the  first  phase,  the  operator  guesses  a  value  for 
nrank  (with  some  guidance  from  Figure  3.4  of  Barber  and  Hill’s 
text),  puts  ntheta  equal  to  twice  nrank,  and  sets  ic*l  to 
instruct  the  program  to  check  for  convergence  over  nrank.  if  the 
solution  converges,  a  message  to  that  effect  is  displayed  on  the 
screen.  Otherwise,  the  process  must  be  repeated  with  another 
value  for  nrank  until  convergence  does  occur.  Once  a 
satisfactory  value  for  nrank  is  established,  a  second  phase  is 
begun  (with  ic=0)  to  check  for  convergence  over  ntheta.  Usually, 
an  ntheta  value  twice  the  size  of  nrank  suffices  to  elicit  the 
"solution  has  converged"  message;  but,  if  not,  larger  values  for 
ntheta  must  be  attempted.  Finally,  in  a  third  session,  the 
operator  enters  the  successful  values  of  nrank  and  ntheta  and 
sets  ic=>2.  The  program  will  determine  a  satisfactory  value  for 
nm  and  then  go  on  to  calculate  and  write  (to  disk)  the  T-matrix. 

Convergence,  and  a  valid  T-matrix,  is  by  no  means 
guaranteed  for  every  set  of  particle  properties  one  might  submit. 
The  program  may  fail  either  because  the  particle  requires  a  value 
of  nrank  >25  (or  ntheta  >100)  or  because  a  run-time  math  error  is 
generated,  which  shuts  the  processing  down.  To  get  a  useful 
gauge  of  when  the  Barber  a*»d  Hill  codes  are  applicable,  we  need 
to  consider  a  goodly  number  of  particles  representing  a  wide 
range  of  sizes,  axial  ratios,  and  refractive  indices,  and  ask  in 
each  case  whether  convergence  can  be  attained  and  with  what 
parameters.  This  could  involve  a  daunting  amount  of  labor  with 
T1  as  it  is  presently  implemented;  so,  we  shall  first  modify  the 
program  to  make  it  run  more  automatically. 


3.  MODIFICATIONS  TO  T1.F0R 

We  wish  to  modify  the  structure  of  T1.F0R,  taking  care 
not  to  alter  its  actual  calculations,  so  that  it  can  be  called  as 
a  subroutine  by  a  new  main  program  we  shall  write  to  take  the 
place  of  the  human  operator.  As  it  comes  from  the  box,  T1.F0R 
comprises  a  main  program  (Tl)  and  seven  subroutines,  most  of 
which  are  called  to  evaluate  special  math  functions.  One 


8 


subroutine,  ADD PRC,  is  of  special  interest  here  because  the 
testing  for  convergence  takes  place  in  it. 

ADDPRC  is  called  twice  when  ic-1  (nrank),  from  lines 
388  and  441,  and  twice  when  ic«0  (ntheta),  both  times  from  line 
388.  (Line  numbers  throughout  this  report  refer  to  positions  in 
the  unaltered  program  T1.F0R.)  In  both  cases,  ADDPRC  writes 
"solution  has  converged"  to  the  screen  and  terminates  the  program 
at  line  1133  -  within  ADDPRC  -  during  its  second  invocation  if 
convergence  has  been  attained;  otherwise,  program  control  returns 
to  Tl,  and  the  program  ultimately  terminates  at  line  453.  The 
screen  message  is  the  only  indication  that  convergence  occurred. 
When  ic-2  (nm),  ADDPRC  may  be  called  numerous  times  from  line 
388,  once  for  each  increment  of  nm,  until  either  convergence  over 
nm  is  attained  or  the  program  crashes  from  a  run-time  math  error. 
Here,  after  writing  the  convergence  message  but  before 
terminating  at  line  1133,  ADDPRC  closes  the  file  't'  (which  was 
written  by  Tl  and  contains  the  now  final  T-matrix)  and  also 
writes  a  file  called  'case'  that  contains  the  particle's  physical 
properties  for  later  identification  by  other  programs. 

To  tell  from  within  the  program  whether  convergence  has 
occurred,  we  defined  a  new  logical  variable  CNOK,  making  it 
common  to  both  Tl  and  ADDPRC.  Just  above  each  line  in  Tl  where 
ADDPRC  is  called,  we  entered  a  new  line  setting  CNOK  false,  and 
at  line  1119  in  ADDPRC  (reached  if  and  only  if  convergence  has 
been  attained)  another  new  line  letting  CNOK  to  true.  We  also 
removed  the  ’stop'  line,  number  1133  in  ADDPRC,  and  put  a 
'return'  at  the  end  of  ADDPRC  so  that  control  returns  to  Tl 
whether  or  not  there  is  convergence.  In  Tl,  following  each  'call 
addprc’  line,  we  now  test  CNOK.  If  it  is  false,  the  program 
keeps  running  as  before;  if  it  is  true,  the  program  jumps  to  the 
'stop'  line  of  Tl  (line  453)  instead  of  being  stopped  from  within 
ADDPRC. 


We  note  here  an  apparent  bug  which,  although  not 
relating  directly  to  the  task  at  hand,  had  to  be  fixed  before 
proceeding.  Lines  1121  through  1126  write  to  the  disk  a  file 
called  'case'  when  the  final  convergence  has  been  achieved,  i.e., 
convergence  over  nm  with  ic«2.  If  a  file  'case'  already  exists 
on  the  disk  -  which  would  be  the  normal  condition,  unless  T1.F0R 
is  being  run  for  the  very  first  time  -  then  the  computer  halts 
operation  with  a  'Stack  Overflow'  error  message  instead  of 
overwriting  the  old  file  with  a  new  'case'  as  intended.  The 
FORTRAN  lines  in  question  appear  to  be  correct,  and  indeed  no 
tinkering  with  them  improved  matters.  Finally,  we  relocated  the 
offending  lines  to  a  new  environment,  in  Tl  after  the  line  (388) 
that  calls  ADDPRC,  since  we  can  now  check  there  for  convergence 
via  CNOK.  The  values  written  into  the  file  'case'  were  already 
known  within  Tl  through  common  blocks  with  ADDPRC.  With  this 
fix,  the  program  wrote  'case'  with  no  further  difficulties.  The 
reason  for  the  problem  was  never  understood  but  is  probably  a 
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quirk  with  the  particular  compiler  used,  Microsoft  Fortran 
Version  5.0. 

The  entire  program  Tl.FOR  was  then  easily  transformed 
to  a  subroutine,  called  T1S,  with  a  few  more  simple  changes.  The 
'stop'  line  453  (which  is  now  always  reached  at  termination)  war 
removed,  and  a  'return'  line  substituted  just  before  the  'end' 
line,  471.  The  read  and  write  lines  for  interfacing  with  the 
operator  were  removed  (lines  76-83  and  105-107),  and  the  whole 
thing  was  named  SUBROUTINE  T1S  (XZ,  AOVRBZ,  CMRZ,  CMIZ,  NRANKZ, 
NTHETZ,  ICZ,  MTXSAX,  CNOKZ,  FIRZ ) .  The  arguments  were  formed  by 
adding  or  substituting  'Z'  for  the  last  letter  of  the  program 
variables  -  new  names  were  required  because  the  variable  names 
already  appear  in  several  common  blocks.  A  new  set  of  statements 
were  written  near  the  beginning  of  T1  to  relate  these  argument 
variables  to  their  originally  defined  corresponding  program 
variables. 


One  last  problem  had  to  be  solved  to  make  this  program 
work  as  a  subroutine.  A  logical  variable  inside  ADDPRC,  called 
FIRST,  is  used  to  tell  whether  ADDPRC  is  being  called  for  the 
first  time  in  any  execution  of  Tl.FOR.  Remember  that  ADDPRC  is 
called  twice  when  ic=l  or  ic*0,  and  numerous  times  when  ic*2,  and 
the  program  branches  differently,  depending  upon  which  invocation 
of  ADDPRC  is  being  processed.  FIRST  is  initially  set  equal  to 
true  at  a  data  statement  (line  926)  each  time  Tl.FOR  is  called  up 
from  the  command  line.  But  when  Tl.FOR  is  made  a  subroutine  to  a 
main  program,  FIRST  will  be  properly  initialized  only  once,  when 
the  whole  program  is  first  started.  All  subsequent  calls  to  the 
subroutine  will  find  FIRST=false,  since  data  statements  are 
executed  only  once,  when  the  program  is  read  in.  To  overcome 
this  problem,  we  put  FIRST  into  a  common  block  shared  by  T1  and 
ADDPRC,  and  included  a  logical  variable  FIRZ  in  the  argument  list 
of  T1S  so  that  FIRST  could  be  properly  set  from  the  main  program 
each  time  T1S  is  called.  A  line  in  T1  sets  FIRST=FIRZ. 

To  control  the  overall  operation,  a  simple  main  program 
called  T1A.F0R  (Figure)  was  written.  The  operator  running  it  is 
asked  to  enter  only  a  particle's  axial  ratio  and  real  and 
imaginary  refractive  indices.  The  program  sets  the  size 
parameter  to  1.0  and  begins  to  check  for  convergence  over  nrank, 
starting  with  nrank*2  and  incrementing  nrank  (called  nrankz  in 
T1A )  until  either  convergence  is  achieved  (indicated  by  CNOKZ 
being  true)  or  until  nrank  is  incremented  above  25,  in  which  case 
the  program  halts.  If  that  convergence  is  successful,  then 
convergence  over  ntheta  is  checked  in  a  similar  fashion,  starting 
with  ntheta  equal  to  twice  nrank  and  increasing,  if  necessary,  in 
steps  of  four.  The  program  halts  if  ntheta  goes  above  100. 
Finally,  convergence  over  nm  is  checked,  and  if  there  is  no 
problem,  the  size  parameter,  nrank,  ntheta,  and  nm  are  written  to 
a  file  on  disk  called  'answer'.  Then,  the  size  parameter  is 
incremented  by  1.0,  and  the  process  begins  again.  This  repeats 
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over  and  over  until  a  size  parameter  is  reached  for  which 
convergence  fails  for  one  of  four  possible  reasons:  nrank 
exceeding  25;  ntheta  exceeding  100;  floating  point  math  error 
M6101 ,  invalid  operation;  or  floating  point  math  error  M6104 
overflow.  These  reasons  for  termination  are  coded  as  R,  T,  1 
and  0,  respectively  in  the  taules  of  results  that  follow.  Once  a 
size  parameter  is  reached  at  which  convergence  fails,  convergence 
will  fail  at  all  larger  size  parameters  (for  the  same  axial  ratio 
and  refractive  index),  though  not  necessarily  for  the  same 
reason . 


A  listing  of  the  program  T1A.F0R  follows: 

PROCRAM  T1A 
LOGICAL  CXOKZ.FIRZ 

C0MM0K  /CMYCOH/  NM , KMV , CMV . TVH , PROOM 
VRITE<«.*>  'EVTER  »o*rb,  «ur,  ■  $ ' 

REAOtS .  •)  A0VR8Z ,  CMRZ  ,  CM  I Z 
0PENIUNIT.8.  FILE-'AKSVER') 

VRITEC8.*) 

VRITEI8,  '  (4I.3A8)  ’  )  ' AOVRB ' , ' MR ' . 'Ml ' 

WRITEC8.  ' (4X ,3F8 .  3) ' )  AOVRRZ.CKRZ.CMIZ 
WHITE! 8.  '  (3X ,4A6) ' )  'X' ,  'mUK'  ,  ’  NTMET  ’  , ' XT  ’ 

DO  28  JM-1.20 
XZ-JRB 


C  . KTJtSAZ  *1  FOR  T'HATMX,  -2  FOR  A- MATRIX 

HTXSAZ  *  1 

C  . . .  CHECK  COHVERGEACE  OVER  HRARX  .  . . 

ICZ*J 

hrarkz-i 


3  KRARKZ  -  KRARKZM 

IFimURKZ.JA.ZS)  THE* 

VRrTE(6,«)  'no  con**Tj*ac«  Icr  nrnak  <  2S' 
COTO  2 
ERD  IF 

RTKETZ-2«RlUmcZ 


FIRZ  •  TRUE. 

CALL  T1S(XZ.AOVR8Z,C»*Z.CXIZ.)OUHK7.RTXETZ. ICZ. HTXSAZ. CHOKZ, FIRZ) 
IF(.NOT.CKOKZ)  COTO  3 

C  .  CHECX  COIfVERGEACE  OVER  HTHETA  . 


ICZ-0 

Kmn-rrHETz-4 

6  RTHETZ*RT KETZM 

irorrum  cT  ioo)  the r 

VRITE16,*)  'no  eonoorjoncn  /or  mthntn  <  $00' 

COTO  21 
END  IF 

FIRZ  -  .TRUE. 

CALL  T1SCIZ,  A0VR8Z  ,CMR2,CJaZ  .RRARKZ,  ITKETZ,  ICZ.MTZSAZ , CRQKZ.FIRZ) 
IF(.HOT.CROia)  COTO  6 


C  .  CHECK  CORYERCERCE  OVER  RM  . 

ICZ-2 

FIRZ  -  .TRUE. 

CALL  TlSUZ.AOTRaZ.CKRZ.OUZ.mm.jrnOTZ.ICZ.fmUAZ.CROKZ.FIRZ) 

C  .  IF  HADE  IT  TO  HERE,  GOOD  CORY,  T-KATRII  WITTER  .... 

WRITE (8, 27)  JM.R7URKZ,  RTKETZ.KXVM 


27  FORMAT  (3X.4I6) 

28  COWTI HUE 
21  STOP 

END 


Figure.  TIA.FOR  PROGRAM 
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4. 


RESULTS 


The  program  T1A  was  run  to  find  convergence  parameters 
for  spheroidal  particles  with  refractive  indices  typical  of 
nonconducting  materials  at  visible  wavelengths  in  air  or  water. 
Tables  1  through  6  show  the  results  for  real  n  equal  to  1.1,  1.4, 
and  1.8,  with  the  imaginary  part  in  each  case  equal  either  to  0.0 
or  to  0.1,  which  is  characteristic  of  a  highly  absorbing 
material.  The  triplet  of  numbers  in  the  clear  boxes  are  the 
convergence  values  of  nrank,  ntheta,  nra  (also  called  nt)  for 
particles  of  the  boxes*  corresponding  size  parameter  and  axial 
ratio.  Shaded  boxes  designate  particles  that  cannot  be  handled 
by  the  Barber /Hill  programs,  and  the  letter  in  the  first  shaded 
box  of  each  column  indicates  the  reason  why  convergence  failed, 
as  described  at  the  end  of  the  last  section. 

Whether  a  solution  converges  for  a  given  size  and  shape 
of  particle  and  the  convergence  parameters  required  to  effect 
that  solution  are  seen  to  depend  only  weakly  on  the  particle's 
refractive  index.  As  either  the  real  or  imaginary  part  of  the 
refractive  index  increases,  the  tendency  is  to  lower  the  size  of 
the  largest  particle,  which  can  be  solved  for  a  given  axial 
ratio,  and  to  increase  the  value  of  nrank  required  to  obtain 
convergence  for  a  given  particle.  For  low  values  of  n,  the 
solvable  particle  with  the  largest  size  parameter  is  a  prolate 
with  an  axial  ratio  of  about  1.5.  With  increasing  n,  the  center 
of  the  distribution  tends  toward  the  sphere,  axial  ratio  1.0. 

The  value  of  nrank  given  in  the  tables  is  the  smallest 
number  that  produces  convergence.  Larger  estimates  may  be  used 
to  run  Tl.FOR;  but,  they  make  the  T-matrix  unnecessarily  large. 

In  fact,  too  large  a  value  of  nrank  may  thwart  convergence.  The 
program  T1A  sets  ntheta  equal  to  twice  nrank  as  an  initial  guess, 
and  in  almost  all  cases,  that  is  adequate.  Often  a  smaller 
ntheta  would  have  sufficed;  but  (unlike  nrank),  there  is  no 
penalty  for  using  a  larger  number.  Values  of  ntheta  even  greater 
than  twice  nrank  are  required  for  some  prolate  particles  with 
very  high  aspect  ratios. 

The  program  was  also  run  with  refractive  indices 
characteristic  of  conductors,  n«0.8  +  il.5,  and  n*1.3  +  i7.1, 
which  represent  gold  and  aluminum,  respectively,  at  a  visible 
wavelength.  The  results,  in  Tables  7  and  8,  are  consistent  with 
the  trends  discussed  above.  In  the  case  of  gold,  we  observed  no 
apparent  problem  with  specifying  a  real  refractive  index  less 
than  1.0. 

One  does  not  have  an  intuitive  feeling  for  the  size  of 
spheroidal  particles  from  their  size  parameters,  and  such  an 
appreciation  would  be  helpful  for  judging  whether  the  Barber/Hill 
programs  are  useful  in  a  given  situation.  Now,  the  size 
parameter  of  a  spheroid  is  defined  to  be  the  product  ka,  where 
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k=2ir/X,  and  a  is  the  rotational  semi-axis.  If  b  is  the 
equatorial  semi-axis,  then  the  volume  of  the  spheroid,  oblate  or 
prolate,  is  given  by  V=(4/3)frab2.  Then,  it  is  easily  shown  that 
the  diameter  of  a  sphere  containing  the  same  amount  of  material 
as  a  spheroid  with  size  parameter  ka  is  given  by 


D  = 


A  ka 
*  (a/b) 2/3 


Using  a  typical  wavelength  for  visible  light,  X  *  0.5  microns,  we 
can  calculate  the  diameter,  in  microns,  of  spheres  equal  in 
volume  to  the  tabulated  spheroids. 

Table  9  presents  the  equivalent  sphere  data.  For 
reference,  a  heavy  line  marks  the  boundary  of  solvable  particles 
with  refractive  index  n-1.4  +  iO.O.  The  size  of  spheroidal 
particles  whose  visible  light  scattering  can  be  calculated  with 
the  Barber/Hill  suite  of  programs  is  surprisingly  small.  At 
best,  spheroids  equivalent  to  spheres  of  about  2  pm  in  diameter 
can  be  treated,  and  for  axial  ratios  much  greater  than  2  (or  less 
than  1/2),  only  particles  equivalent  to  submicron  spheres  can  be 
handled . 


5 .  CONCLUSIONS 

The  information  generated  and  displayed  in  Tables  1 
through  8  should  help  users  of  Barber  and  Hill’s  T-matrix  codes. 
The  tables  indicate  the  size  and  shape  range  of  those  spheroidal 
particles  whose  light  scattering  can  be  calculated  with  the 
codes,  and  specify  trial  values  for  convergence  parameters  so 
that  Tl.FOR  can  be  run  quickly  and  efficiently. 

To  save  time  running  the  many  cases  needed  to  produce 
the  tabulated  data,  we  wrote  a  modified  version  of  Tl.FOR  which 
can  accept  particle  properties  as  input  and  write  the  particle's 
T-matrix  as  output  with  no  further  intervention  by  a  human 
operator,  as  required  with  the  original  Tl.FOR.  That  suggests 
another  variation  on  the  Barber  and  Hill  codes.  If  an  automated 
T-matrix  writing  routine  can  be  integrated  with  a  basic  light 
scattering  calculation  code  such  as  T3.F0R,  then  it  should  be 
possible  to  produce  an  independent  subroutine  which,  when  given 
particle  properties  and  orientation,  returns  a  fundamental 
scattering  property.  A  small  supply  of  such  subroutines  would 
allow  a  casual  programmer  to  write  easily  customized  programs  to 
supplement  the  applications  already  provided  by  Barber  and  Hill. 
We  plan  to  explore  this  possibility  in  the  near  future. 
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TABLE  1 .  Convergence  values  of  (nrank,  ntheta,  nt)  for  n=l.  1  +  *0.0 
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TABLE  2.  Convergence  values  of  (nrank,  ntheta,  nt)  for  n=  1.1  +  iO.l 
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TABLE  3.  Convergence  values  of  (nrank,  ntheta,  nt)  for  n=1.4  +  iO.O 
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TABLE  4.  Convergence  values  of  (nrank,  ntheta,  nt)  for  n«1.4  +  iO.l 
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TABLE  5.  Convergence  values  of  (nrank,  ntheta,  nt)  for  n»1.8  +  iO.O 
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TABLE  6.  Convergence  values  of  (nrank,  ntheta,  nt)  f  x  n«1.8  +  iO.l 
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TABLE  7.  Convergence  values  of  (nrank,  ntheta,  nt)  for  n=0.8  +  ii.5 
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